####################
# Replication codes for "Does Exposure to Election Fraud Research Undermine Conﬁdence in Elections?"
# 2. Figures and Table
# Last modified: 2024-03-10
####################

### install packages if not already installed
# install.packages("here")
# install.packages("broom")
# install.packages("tidyverse")
# install.packages("stargazer")

### set working directory here
# setwd("")

library(here)
library(broom)
library(tidyverse)
library(stargazer)

df <- read_csv(paste(here(),"data/election_fraud_survey_cleaned.csv", sep="/"))

# select variables for descriptive plot
# create labels for descriptive plot variables
df_desc <- df %>% filter(mcheck == 1) %>% select(core_dv2, dv1, dv2, clicked) %>% 
  mutate(clicked = case_when(clicked == 0 ~ "Clicked Nothing",
                             clicked == 1 ~ "Demand Investigation", 
                             clicked == 2 ~ "Oppose Investigation"),
         dv1 = ifelse(dv1 == 1, "Yes", "No"),
         dv2 = ifelse(dv2 == 1, "Yes", "No"),
         core_dv2 = ifelse(core_dv2 == 1, "Yes", "No"))

# create a long format data frame
df_desc <- df_desc %>% pivot_longer(everything(), names_to = "vars")

# describe variable names
# reorder factors
df_desc <- df_desc %>% mutate(vars = case_when(vars == "core_dv2" ~ "Do you think fraud occurred in the election?",
                                     vars == "dv1" ~ "Do you trust the election results?",
                                     vars == "dv2" ~ "Do you think the election results reflected the will of the people?",
                                     vars == "clicked" ~ "Petition Link Clicked")) %>%
  drop_na(value) %>%
  mutate(vars = fct_relevel(vars, "Do you think fraud occurred in the election?", "Petition Link Clicked"),
         value = fct_relevel(value, "Yes", "No", "Clicked Nothing", "Demand Investigation", "Oppose Investigation"))

#### Figure 1
pdf("figure/fig1_descriptive.pdf", 10, 6)
df_desc %>% count(vars, value) %>% group_by(vars) %>% mutate(outcome = n/sum(n) * 100) %>%
  ggplot(aes(value, outcome, fill = value)) +
  geom_col(color = "gray") + facet_wrap(~ vars, scales = "free_x") +
  theme_bw() +
  scale_fill_manual(values = c("#386cb0", "#ffff99", "#f1a340", "#f7f7f7", "#998ec3")) + 
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(size = 11),
        strip.text = element_text(size = 12)) +
  labs(x = "",
       y = "Percent")
dev.off()

#### Figure 2: Fraud Belief

# run regression 
m1 <- lm(core_dv2 ~ treat_binary, df %>% filter(mcheck == 1))

# create data frame for figure 2
m1 <- m1 %>% tidy() %>%
  mutate(ul = estimate + 1.96 * std.error,
         ll = estimate - 1.96 * std.error) %>%
  slice(-1) %>%
  mutate(term = gsub("treat_binary", "", term),
         term = factor(term, levels = c("Unlikely",
                                        "Likely")),
         term = recode(term, 
                       "Likely" = "Substantial Chance",
                       "Unlikely" = "Slim Chance"))

# create figure 2
pdf("figure/fig2_fraud_belief.pdf", 4, 4)
ggplot(m1, aes(x=term, y=estimate, ymin=ll, ymax=ul)) +
  geom_pointrange() + theme_bw()  +
  coord_flip() + xlab('') + ylab('Treatment Effect') +
  #facet_wrap(~ type, ncol=1, scales = "free_y", strip.position="right") + 
  geom_hline(aes(yintercept=0), linetype = "dashed", color = "gray") +
  theme(axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        strip.text = element_text(size = 12)) + ggtitle("Belief about Election Fraud") +
  ylim(-0.3, 0.3) 
dev.off()






#### Figure 3: Fraud Belief by PID

# run regression
model1 <- lm(core_dv2 ~ moon_bin * treat_binary, df %>% filter(mcheck == 1))

# create data frame for figure 3
m1 <- model1 %>% tidy() %>% slice(-1:-2) %>%
  mutate(term = gsub("treat_binary", "", term),
         pid = ifelse(grepl("moon_bin", term), "Moon Supporter", "Moon Non-supporter"),
         term = gsub("moon_bin:", "", term)) 
m1$estimate[3] <- m1$estimate[1] + m1$estimate[3]
m1$estimate[4] <- m1$estimate[2] + m1$estimate[4]
m1$std.error[3] <- sqrt(vcov(model1)[3,3] + vcov(model1)[5,5] + vcov(model1)[3,5])
m1$std.error[4] <- sqrt(vcov(model1)[4,4] + vcov(model1)[6,6] + vcov(model1)[4,6])
m1 <- m1 %>% 
  mutate(ul = estimate + 1.96 * std.error,
         ll = estimate - 1.96 * std.error,
         likely = factor(term, levels = c("Unlikely", "Likely")),
         likely = recode(likely,
                         "Likely" = "Substantial Chance",
                         "Unlikely" = "Slim Chance"),
         pid = fct_recode(pid, "Approve of President" = "Moon Supporter", 
                          "Disapprove of President" = "Moon Non-supporter"))


# create figure 3
pdf("figure/fig3_fraud_belief_partyid.pdf", 9, 7)
ggplot(m1, aes(x=likely, y=estimate, ymin=ll, ymax=ul, color = pid, shape = pid)) +
  geom_pointrange(position = position_dodge(width = -0.5), size = 0.7, lwd = 0.7) + theme_bw()  +
  coord_flip() + xlab('') + ylab('Treatment Effect') +
  #facet_wrap(~ type, ncol=1, scales = "free_y", strip.position="right") + 
  geom_hline(aes(yintercept=0), linetype = "dashed", color = "gray") +
  theme(axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 13),
        strip.text = element_text(size = 13)) + ggtitle("Belief about Election Fraud by Party ID") +
  ylim(-0.35, 0.35) + facet_wrap(~ pid) +
  labs(color = "Presidential Approval",
       shape = "Presidential Approval") +
  scale_color_manual(values = c("#d7191c", "#2b83ba")) +
  scale_shape_manual(values = c(19, 15)) +
  theme(legend.position="right")
dev.off()


#### Figure 4: behavioral measure

# DVs for behavioral measures
df <- df %>% mutate(click_demfraud = ifelse(clicked == 1, 1, 0),
                    click_futurefraud = ifelse(clicked == 2, 1, 0))

# run regression for each behavioral measure
m1 <- lm(click_demfraud ~ treat_binary, df %>% filter(mcheck == 1))
m2 <- lm(click_futurefraud ~ treat_binary, df%>% filter(mcheck == 1))

# create data frame for figure 4
m1 <- m1 %>% tidy() %>%
  mutate(ul = estimate + 1.96 * std.error,
         ll = estimate - 1.96 * std.error) %>%
  slice(-1) %>%
  mutate(term = gsub("treat_binary", "", term),
         term = factor(term, levels = c("Unlikely", "Likely")),
         type = "Demand Investigation",
         likely = c("Likely", "Unlikely"))


m2 <- m2 %>% tidy() %>%
  mutate(ul = estimate + 1.96 * std.error,
         ll = estimate - 1.96 * std.error) %>%
  slice(-1) %>%
  mutate(term = gsub("treat_binary", "", term),
         term = factor(term, levels = c("Unlikely", "Likely")),
         type = "Oppose Investigation",
         likely = c("Likely", "Unlikely")) %>%
  bind_rows(m1) %>%
  mutate(likely = recode(likely,
                         "Likely" = "Substantial Chance",
                         "Unlikely" = "Slim Chance"))



# create figure 4
pdf("figure/fig4_behavioral.pdf", 6, 5)
ggplot(m2, aes(x=likely, y=estimate, ymin=ll, ymax=ul)) +
  geom_pointrange() + theme_bw()  +
  coord_flip() + xlab('') + ylab('Treatment Effect') +
  #facet_wrap(~ type, ncol=1, scales = "free_y", strip.position="right") + 
  geom_hline(aes(yintercept=0), linetype = "dashed", color = "gray") +
  theme(axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 13),
        strip.text = element_text(size = 13)) + ggtitle("Behavioral Measures: Click Petition Links") +
  facet_grid(~ type) +
  ylim(-0.3, 0.3)
dev.off()


#### Figure 5: bush_prather DVs

# run regression for each bush and prather DV
m1 <- lm(dv1 ~ treat_binary, df %>% filter(mcheck == 1))
m2 <- lm(dv2 ~ treat_binary, df%>% filter(mcheck == 1))

# create data frame for figure 5
m1 <- m1 %>% tidy() %>%
  mutate(ul = estimate + 1.96 * std.error,
         ll = estimate - 1.96 * std.error) %>%
  slice(-1) %>%
  mutate(term = gsub("treat_binary", "", term),
         term = factor(term, levels = c("Unlikely", "Likely")),
         type = "Confidence",
         likely = c("Likely", "Unlikely"))


m2 <- m2 %>% tidy() %>%
  mutate(ul = estimate + 1.96 * std.error,
         ll = estimate - 1.96 * std.error) %>%
  slice(-1) %>%
  mutate(term = gsub("treat_binary", "", term),
         term = factor(term, levels = c("Unlikely", "Likely")),
         type = "Reflect Will of People",
         likely = c("Likely", "Unlikely")) %>%
  bind_rows(m1) %>%
  mutate(likely = recode(likely,
                         "Likely" = "Substantial Chance",
                         "Unlikely" = "Slim Chance"))



# create figure 5
pdf("figure/fig5_ther.pdf", 6, 5)
ggplot(m2, aes(x=likely, y=estimate, ymin=ll, ymax=ul)) +
  geom_pointrange() + theme_bw()  +
  coord_flip() + xlab('') + ylab('Treatment Effect') +
  #facet_wrap(~ type, ncol=1, scales = "free_y", strip.position="right") + 
  geom_hline(aes(yintercept=0), linetype = "dashed", color = "gray") +
  theme(axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 13),
        axis.title.x = element_text(size = 13),
        strip.text = element_text(size = 13)) + ggtitle("") +
  facet_grid(~ type) +
  ylim(-0.3, 0.3)
dev.off()


#### Table 1


# run regression models 1-6
t1 <- lm(core_dv2 ~ treat_binary, df %>% filter(mcheck == 1))
t2 <- lm(click_demfraud ~ treat_binary, df %>% filter(mcheck == 1))
t3 <- lm(click_futurefraud ~ treat_binary, df%>% filter(mcheck == 1))
t4 <- lm(core_dv2 ~ moon_bin * treat_binary, df %>% filter(mcheck == 1))
t5 <- lm(dv1 ~ treat_binary, df %>% filter(mcheck == 1))
t6 <- lm(dv2 ~ treat_binary, df%>% filter(mcheck == 1))

# create table 1
stargazer(t1, t4, t2, t3, t5, t6, report=('vc*p'),
          dep.var.labels = c("Belief about Fraud",
                             "Demand Investigation",
                             "Oppose Investigation",
                             "Confidence",
                             "Reflect Will"),
          order = c(2,3,1,4,5,6),
          covariate.labels = c("Substantial Chance",
                               "Slim Chance",
                               "Approve Moon",
                               "Approve $\\times$ Substantial",
                               "Approve $\\times$ Slim",
                               "Constant"),
          no.space = TRUE, 
          keep.stat = c("n", "rsq", "adj.rsq"),
          star.cutoffs = NA,
          title = "Treatment Effects Across Outcome Variables")

